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Abstract 

Using the transfer matrix method, we give the exact solution of a determin- 
istic sandpile model for arbitrary N, where is the size of a single toppling. 
The one- and two-point functions are given in term of the eigenvalues of an 
N X N transfer matrix. All the n-point functions can be found in the same 
way. Application of this method to a more general class of models is discussed. 
We also present a quantitative description of the limit cycle (attractor) as a 
multifractal. 
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I. INTRODUCTION 



Sandpile model has been conceived as the simplest model which can illlustrate self- 
organized criticality(SOC). It has developed into an interesting cellular automata model 
of nonlinear dynamics in its own right. Bak, Tang, and Wiensenfeld are the first to 
use numerical method to study some of the sandpile models and observed that the models 
automatically evolve into a self-organized critical state while they possess 1// spectra in 
both spatial and temporal distributions of certain physical quantities. It was conjectured 
that SOC maybe an universal characteristic underlying the nonlinear, dispersive systems in 
nature such as earthquakes, forest fire, turbulence, etc... 

The simplest cellular automata model of sandpile is to assign a height number hi to each 
site on a one dimensional lattice with length L. There are two basic operations of the model 
- dropping and toppling. Dropping means that one sand is added at some site of the lattice, 
i.e. Di : hi ^ hi + 1. Toppling occurs when a slope (defined as the difference in height 
between adjacent sites) exceeds some critical value, A^. If toppling occurs at one site, some 
sands at the site will be moved to the other sites which may trigger further topplings. 

Even with these simple rules of evolution, the analytic solution of the sandpile model are 
typically very difficult to come by especially when L or N becomes large, or when the lattice 
dimension is more than one. Most researcher handle the models by numerical simulation 
P-^. If the rules are such that the evolution of the system is independent of the order of the 
droppings, then the model is called abelian. For a large class of abelian models some exact 
results have been obtained by Dhar et al. ^ The non-abelian models are typically much 
harder to solve. Here we shall investigate a class of nonabelian models in the deterministic 
limit, i.e., the sand being dropped at a fixed site. 

We consider the one dimensional case and label the sites from left to right as 1 to L. The 
sand is dropped only at the site 1. If the slope at a site exceeds a given number A^, then 
the sand will topple to the right. Let the slope ai = hi — /ij+i, then the toppling rule is "if 
(Tj > N, then (Tj_i a^.i + N , ai ^ ai — {N + 1) , and ai+N — > ai^N + 1" • The rule should 
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be modified when toppling occurs near the boundary. The condition at the left boundary 
is trivial. When sands reach beyond the right boundary they drop out from the system, i.e. 
we keep hi = for i > L. A state in which all aj < is called a stable state. Toppling stops 
when a stable state is reached. Each dropping and subsequent toppling processes will result 
in the transition from one stable state to another. Since there are only a finite number of 
different states in the system, after dropping enough sand at site 1, the system will step into 
a cycle called the limit cycle. For the system that we consider here there is only one limit 
cycle in the problem. Not all the stable states are in the limit cycle. Those in the limit cycle 
are called allowed states. In the following, all the states we will refer to are allowed states. 
The number of allowed states in the limit cycle is A^^. 

Some special solvable non-abelian models have been investigated in Ref. [0-0. In Ref. 
0, the simplest (A^ = 2) of a class of non-abelian sandpile models was solved. Some general 
property of the model for arbitrary was also exposed including some results about the 
random drop cases |p. In Ref. a new method using the transfer matrix idea from the 
statistical mechanics models was developed and applied to solve the A^ = 3 case of the same 
class of sandpile models. It was commented that the method should also be useful for solving 
the models with arbitrary A^. In this note, we wish to develop this idea further and supply 
a solution of the model with arbitrary A^. 

The method to be introduced here is particularly useful for those models in which the 
structure of the limit cycles has been worked out. Once the cycle structure is known, its 
information can be succinctly summarized in a matrix which we call the transfer matrix 
by analogy with the similar matrix in statistical mechanics. For the deterministic model 
defined above, the structure of the limit cycle has been worked out in ref. [0]. Therefore we 
should use it as our main example, the method may be applicable to a much wider classes 
of models. A multifractal description of the limit cycle is also worked out. 
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II. THE DEFINITION OF TRANSFER MATRIX 

A allowed states of length L in the hmit cycle are characterized by three conditions: 

1. (7i > for < i < L. 

2. There exists at least one site i for any consecutive N sites such that Ui — A?", (stability 
condition) 

3. There exists a site i satisfying L — ul < i < L such that (7j = A^. (boundary condition) 

If [(Ji(T2...(7l] is an allowed state in the limit cycle with length L, then [(Ti+i(7i+2---CL] is 
also an allowed state in the limit cycle with length L — i. 

A. Reduced basis 

To make the method of tranfer matrix more powerful, it is useful to define a reduced 
basis, I \n)i, which classifies the allowed states according to the constraint on their preceding 
sites (to their left). The set of allowed states of length / can be classified according to a 
quantum number n, 1 < n < N. A state [r)irj2---r)i\, which is an allowed state in the limit 
cycle with length is said to have the quantum number n if rjj^+i^n — N and rji < N for 
i < N + 1 — n when A^ + 1 — n < /, while rji = I + n — 1 and r]i < N ior i < I when 
N + l—n > I. A collection of such states can be denoted as \ \n)i and called a reduced state, 
while \n)i denotes a typical state in ||n);. Every allowed state belongs to one of the reduced 
states. We therefore call the set of \\n)i, with 1 < n < A^, a reduced basis. The quantum 
number n labels the fact that, for any state [r)irj2---rji\ in one of the n slopes to the left 
of r]i has to be A". 

For a allowed state of length L, [aia2...<JL], the last / slopes [(Ti,_;+i(TL_/+2...cri] is an 
allowed state of length I, which belongs to some reduced basis \n)i, and [(Tia2---(TL-icrL^i+i{— 
n)] is an allowed state of length L — l + 1. This correspondence is one-to-one for fixed I and 
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this provides a route to construct the hmit cycle of length L from the hmit cycle of length 
l{< L). 



B. The time-average of space dependent functions 

Let Si be the number of allowed states of length I + 1 with (T;+i = N. Prom the three 
conditions for allowed states in the limit cycle, we have Si = {N + 1)' for < / < A^, 
Sj^ = {N + 1)^ - N^, and Si+i = {N + l)Si - N^Si_n for / > N. For any allowed state of 
length L with (jj = N, [(7j+i(Tj+2---CL] can be any allowed state of length L — i. The number 
of allowed states of length L with ai = N is thus N^~^Si^i. 

It should be noted that the three conditions for allowed states in the limit cycle only 
differentiate ai — N and N unless i — L. Since we know the total degeneracy of states 
for (Tj taking arbitray value and the degeneracy for states in which = N can be calculated 
using Si, the degeneracy for ^ N can be easily obtained. Therefore, Si contains all 
the information about the time average of space dependent functions, such as the n-point 
function < ai...aj >. Here < ... > means the average in time over a limit cycle. The one- 
point function can be evaluated in terms of the summation of all allowed states in a limit 
cycle, i.e., 

NN^-'Si-i + ^(N^ - N^-^Si_i) 
= k 

= V + ^<^- (2-1) 

Similarly, the two-point function is 

(TV -1)2 N^-l fSi^, , Vi 



< (T,(T-,- >= 1 : h 

' ^ 4 4 V A^' A^^ 
(N + irS^S^ 
+ 4 W ' ^<2<1^- (2.2) 

Actually, all the n-point functions can be expressed in a similar form. Therefore, the calcu- 
lations of the n-point functions are reduced to the calculation of Si. It should be noted that 
these time-average of space dependent functions are independent of L. 



One can use mathematical induction to get the following close form for Si. 



SmN+l={N+l) 



mN+l 



\m - j)N + 1^^^^^^^ ({m - j)N + 



(2.3) 



with m > and < I < N. However, such tedious close form is difficult to use. A better 
form can be obtained by transfer matrix method. 

C. Transfer matrix 



For any state in \n)i, we can add a preceding slope a — N, resulting in a state in 
\N)i+i, or add a preceding slope a = 0,1,2, N — 1 if n > 1, resulting in a state in 
\n — We can express this property in term of a transfer operator T such that T| \n)i = 

||A'");+i + N\\n — if n > 1, and T\\l)i — ||A^);+i. For example, the first equation 

represents the fact that, given a typical state, [riiri2...rii], in \ \n)i the corresponding allowed 
state [f], r]ir]2...r)i] contains one state in ||A'")i_|_i and N states in \ \n — As far as transfer 

matrix is concerned, all the allowed states within the same reduced state are identical and 
there is no need to distinguish them. It is also clear from the definition of transfer matrix 
that its properties are independent of the subscript /. Since it is obvious that each operation 
of T increases the I by one, the subscript can be ignored all together. 

Represent \\N) by col{l ... 0), \\N - 1) by col{0 1 ... 0), ||1) by col{0 ... 1), the 

transfer operator can thus be represented as an x matrix: 

/I 1 1 ... 1 1\ 

... 

AT ... 





V 





A^ o; 



(2.4) 



Note that co/(l 1 .. 
length one are ai - 



. 1) is an eigenvector of T of eigenvalue A^. The allowed states of 
= 1,2,..., A?", which correspond to the reduced basis ||l)i, ||2)i, 
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||A^)i, respectively. The number of allowed states of length L is easily proved to be 
(1 1 ... l)T^-^col{l 1 ... 1) = A^-^. This is the simplest way of deriving this result pointed 
out already in ref. |^. The number of states corresponds to the reduced basis \n)i with 
1 < n < is {5N,n 5N-i,n ■■■ 5i,n)T^~'^col{l 1 ... 1) = N^~^ , which is independent of ra, and 

Si = {ll ... 1)TW(1 ... 0). (2.5) 

Si can now be expressed in terms of the eigenvalues of the transfer matrix T, and we 
have 

yyZ+l / N-1 \ 

where A^Ao(Ao = 1), A^Ai, NX^^i are the eigenvalues of T and Ai, A2, Aat-i are the 
solutions of the equation 

1 + 2A + SA^ + ■ • • + NX^'^ = 0. (2.7) 

The one and two-point functions are thus 

N+1 1^-^ 

< >= - 2 E (2.8) 

m=l 

and 

< >= . + E + E 

^ ^ m=l ^ m=l 

■1 Af-1 1 N-l N-l 

+^ E + 7 E ^™ E i<3<L. (2.9) 

^ m=l ^ m=l k=l 

For N = 2 and 3, the results agree with [^0. 

With y = 1/X and multiplying Eq. (pl7| ) by ?/ — 1, we have 

^ ^AT-l + . . . + ^ = ]v ^ |y| > 1 ^ |A| < 1 (2.10) 
where the redundant solution y = 1 is excluded. Multiplying hj y — 1 again, we have 

yN+i = + _ ^ + 1)1^1 + > |^|iv+i^ (2.11) 



which imphes \y\ < 2 ioi N > 2 and therefore \y\ < {3N + 2)^/'^^"'"^). The solutions of y 
are inside the annulus 1 < \y\ < ymax = (3A^ + 2)^/*^^+^^ In fact for large one can show 
that ymax = 1 + N^^log{2N) and therefore as ^ oo all the roots have \y\ near 1. From 
Eq. ( p.ll|) we have 

N 



yr+^ = (iv+i) 



y 



7V+ 1 ' 

which means the solutions with smaller \y\ are all near y = 1 for large A^. 

To explore the asymptotic behavior of 5*;, we are interested in the solutions of Eq. ( p.7| ) 
with larger |A|. Therefore, we would like to find the solutions of Eq. with smaller 

\y\ excluding the redundant solution y = 1. Considering only one of the complex conjugate 
solutions in the upper half of the complex plane, let yj = a^e*^^ with < j < A^/2 and 
(f)j = i[(N + l)yj - N] with < (f)j < re such that (A^ + 1)^^- = + 2j7r. It should be noted 
that (A^ + 1)6*0 = 00 results in the redundant solution yo = 1. From Eq. (|2.11| ) we have 



^2N+2 = (jv + l)2a| + N^- 2N{N + l)a,- cos^, 

and 

(A^ + l)a,- sin6',- 

tan<7)i = - — ^ ^ — — — 

{N + l)ajCosej - N 

For large A^ we only have to consider the case of j/N <^ 1 to find the smaller \y\ solutions. 
Let aj = 1 + Xj/{N + 1) and Oj = (2j7r + (f)j)/{N + 1) < 1, we have (as N ^ oo) 

e2-. _ (1 + a;.)2 ^ (0^. + 2jnf (2.12) 

and 

2j7r + (hi 

tsin(j)j = — (2.13) 

1 + Xj 

where the trigonometric functions have been expanded to the second order. Table | shows 
the values of the first few x's and 0's as A^ — > oo. From Eq. ( [^.9D , the correlation length is 
found to be N/xi for large A^. 
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III. MULTIFRACTAL CHARACTERIZATION OF THE LIMIT CYCLE 



Now we would like to describe the limit cycle (attractor) quantitatively. In general, the 
state in the limit cycle corresponds to a path in a tree. It is now possible to identify a 
multifractal along the time axis. Rescaling the time variable so that each time step is of 
length S — N"^, the duration for a limit cycle fits into the time interval [0, 1]. Consider 
the states [aia2...(JL\ in the limit cycle which contains Nq O's, Ni I's, and Ntv A'"'s,. 
The corresponding time steps for such states form a fractal in the interval [0, 1] with fractal 
dimension /(^o, Ci; Cat) where = Nj/L, for < i < and satisfy Y.i — 1- The above 
relations are exact as L ^ oo. 

For N — 2/\i is straight forward to get 



/(Co, a, 6) ^ 



In 2 



Coin i ^'-^r^' \ +61n f ^'-^r^' \ +61n ^ 



(3.1) 

This function depends on and ^iva. the same way, which suggests that the boundary 
effect due to the specific constraint on ctl is negligible in the large L limit. Neglecting this 
boundary effect, there is no more difference between n — 0, 1, and N — 1. Thus it is 
natural to define the fractal dimension as a function of ^ — Nf^/L with < ^ < 1. The 
problem is equivalent to counting the ways to distribute (1 — objects into slots and 
every slot admits at most N — 1 objects. Now we have 



P^l-Nn (3.2) 

x=0 



(L-NAr)!dx^-N^ \x-l J 
The is as L — > oo. After taking the differentiation with x apart on [(x^ — l)/{x — 



N-l 



jyfH^-Of ^ J2 Ar'=(i-/)+*^^/' (3.3) 

fe=0 

where /(^ + is approximated as /(O + A^/'. Let f{^)=g{^)/lnN+l-^ with R = e^' 
and Z = e~^^~^^'\ which means 



(3.4) 



Now Eq. ( ^.3[ ) becomes 

R = Z + Z'^ + --- + Z^. (3.5) 
It follows from the definitions of R and Z that 

g{i) = i HZ + Z2 + . . . + Z^) - In Z (3.6) 

with 



Z + Z"' 



7N-1 



^ 1 + 2Z + 3Z^ + --- + NZ^-^' ^^■''^ 

This result obviously satisfy the boundary conditions /(^ = 1/A^) = 1 — and /(^ = 
1) = 0. /(^) = 1 and /'(O = happen at ^ = 2/{N + 1), which is consistent with 



< ai >^ {N + 1) /2 as i ^ oo from Eq. (|2.8|) . Equations ( |3.6|) and (|3.7| ) express the fractal 
dimension /(^) in terms of the variable Z , < Z < oo. 



IV. DISCUSSIONS 

We have introduced the concept of transfer matrix into the study of steady properties 
of sandpile models. This is closely related to the Hamiltonian formulation of the usual 
statistical mechanics. Indeed, one may regard the formulas for the correlation function such 
as equation ( |2.1| ) and ( p.2|) as a "path integral" expressions in a discrete formulation. The 
transfer matrix plays the role of the evolution operator. 

The usefulness of the transfer matrix formulation was illustrated by deriving the one- and 
two-point correlation functions for a deterministic sandpile model with an arbitrary critical 
slope N. We found that the two-point function decreases exponentially as the separation 
of the two point increases with a correlation length depending on A^. In the unit of lattice 
spacing, the correlation length is found to be the largest (in magnitude) root of an (A^-l)th 
order polynomial. The correlation length is proportional to in the large limit. 
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In the L — > oo limit, the states in the hmit cycle can be characterized as a multifractal. 
We derive an explicit expression for the multifractal dimension for arbitrary A'^. 

Even though there exists only short range spatial correlation for any N, power law 
behavior is expected for the time-time correlation function. Derivation of an explicit formula 
for this correlation function is under investigation. 
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TABLES 

TABLE L Somo numerical values of x . and r>, as — > oc. 



i 

1 2.089 1.178 

2 2.664 1.312 

3 3.026 1.374 

4 3.292 1.410 

5 3.501 1.435 
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